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We study the low-energy collective oscillations of a dilute Bose gas at finite temperature in the 
collisionless regime. By using a time- dependent mean-field scheme we derive for the dynamics of the 
condensate and noncondensate components a set of coupled equations, which we solve perturbatively 
to second order in the interaction coupling constant. This approach is equivalent to the finite- 
' temperature extension of the Beliaev approximation and includes corrections to the Gross- Pitaevskii 

I theory due both to quantum and thermal fluctuations. For a homogeneous system we explicitly 

^\ ■ calculate the temperature dependence of the velocity of propagation and damping rate of zero sound. 

In the case of harmonically trapped systems in the thermodynamic limit, we calculate, as a function 
^ , of temperature, the frequency shift of the low-energy compressional and surface modes. 
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, The study of collective excitations is one of the main areas of interest for the experimental and theoretical research 
J — ■ activity in trapped Bose-condcnsed gases (for a review of experimental and theoretical investigations see respectively 
^f) and 1^). At low temperatures, the frequencies of the low-energy collective oscillations of the condensate have been 
T-H ' measured with great accuracy and found in very good agreement with the predictions of the mean-field Gross- 
[ Pitaevskii theory |^,|| . In a series of experiments carried out at JILA [Q and MIT |^ the excitations of a trapped Bose 
■ gas have also been explored as a function of temperature. The main features are: on the one hand oscillations of both 
the condensate and the thermal cloud are visible and, on the other hand, the oscillations are increasingly damped as 
temperature is raised and temperature dependent frequency shifts are also observed. A theoretical description which 
correctly accounts for these phenomena has not yet been fully developed. 
5^ , At finite temperature the dynamics of Bose-condensed systems is complicated. Depending on the temperature, 
= density and frequency of the excitations one is probing different regimes (for an exhaustive discussion see the books 
^ [ 1^ and ||lo[| ). If the frequency uj is much smaller than the inverse of the typical collision time Tc. luTc ^ 1, the 
O ■ excitations are collective collisional modes, which are described by the theory of two-fluid hydrodynamics. In terms 
^ ' of length scales this regime is equivalently defined by the condition: X^x S> £mfp, where Xex is the wavelength of the 
^ , excitation and imfp is its mean free path. At low temperatures and low values of the density the mean free path 
becomes comparable with the size of the system. In this case, which corresponds to the condition loTc 3> 1, one 
^ . is probing the collisionless regime, which is properly described by mean-field theories. Collisionless modes can be 
further distinguished into collective and single-particle excitations, depending on whether the excitation energy lies 
respectively well below or above the chemical potential /i. Single-particle excitations have wavelength much smaller 
than the healing length of the condensate, which is defined as ^ = where a is the s-wave scattering length 

and no is the condensate density. On the contrary, collective modes satisfy the condition: X^x ^ C- Finally, in 
harmonically trapped systems, collective modes can behave semiclassically if their energy is much larger than the 
typical trapping energy: hujho ^ <C fi, where tOho is the harmonic oscillator frequency. If instead Hlu ^ hcuho, the 
discretization of levels becomes important and one is not allowed to treat the excitation as quasiclassical. 

Collective modes in the collision-dominated regime have been investigated in harmonically trapped systems by 
several authors . The present work is focused on the study of collective excitations in the collisionless regime. In 
the last years a large number of theoretical papers have appeared in the literature addressing this problem. Mean-field 
approaches, which extend to finite temperature the Gross-Pitaevskii equation for the order parameter, have been put 
forward and applied to the calculation of the low-energy modes in traps ]T3|-p^. However, in these approaches 
the noncondensate component is treated as a static thermal bath and its dynamic coupling to the oscillations of 
the condensate is neglected. The results obtained for the collective modes do not adequately reproduce the features 
observed in experiments, in particular these approaches do not account for the damping of the excitations. 

More accurate time-dependent mean- field schemes have been proposed ]l^[l9t , which describe the coupled dynamics 
of the condensate and noncondensate components. These methods have been applied to the study of damping in 
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trapped systems and agree with results obtained from perturbation theory jEQ£l| . Exphcit calculation of the 
damping rate of the low-energy modes in harmonic traps has been carried out in [pl|-|24[ and found in good agreement 
with experiments. Similar methods have also been applied to the calculation of the temperature dependence of the 
frequency shifts in the coUisionless regime. Bijlsma and Stoof have used a variational ansatz to describe the 
time evolution of the condensate and the thermal cloud and have calculated the frequencies of the coupled modes in 
which the two components move either in phase or out of phase. These authors also suggest that the avoided crossing 
between the in and out of phase modes might be the reason of the features observed for the m = mode at JILA . 
Olshanii has explicitly analyzed the JILA m — 2 mode 0, suggesting that the observed temperature dependence of 
the excitation frequency might be due to a strong resonance between the oscillation frequency of the condensate and 
one of the eigenfrequencies of the thermal cloud . Fedichev and Shlyapnikov ||2^ have developed a Green's function 
perturbation scheme for inhomogeneous Bose-condensed gases at finite temperature and have calculated energy shifts 
and damping rates of quasiclassical collective modes, which satisfy the condition TiLDho e fJ., being e the energy 
of the excitation. Very recently, Reidl et al. have calculated by the dielectric formalism the frequency shift of the 
m — and m — 2 modes, and compared the results with the JILA experiments. 

In the present work we derive, within a time-dependent mean-field scheme, coupled equations for the dynamics 
of the condensate and noncondensate components. These equations are solved perturbatively to second order in 
the interaction coupling constant. For homogeneous systems this approach is equivalent to the finite-temperature 
extension of the Beliaev approximation discussed in psf . In the homogeneous case we give explicit results for the 
temperature dependence of the velocity of zero sound, which include effects beyond the Bogoliubov theory. We also 
apply our analysis to harmonically trapped systems in the thermodynamic limit. In this regime, which is reached 
for systems with a very large number of trapped particles, one can use the Thomas-Fermi approximation for the 
condensate and neglect finite-size effects. Under these conditions, which are not difficult to realize in experiments 
(see e.g. |^), the frequencies of the collective modes are found to change with temperature due to static and dynamic 
correlations beyond the Gross-Pitaevskii theory. We calculate, as a function of temperature, the frequency shifts of 
the lowest compressional and surface modes. We find that at the intermediate temperatures T ~ 0.6-0.7 Tc, where Tc 
is the transition temperature, the fractional shift due to beyond mean-field effects is of the order of 5%. This result 
should be compared with the corresponding correction predicted at very low temperatures |^,^ and arising from 
quantum fluctuations, which turns out to be typically of the order of 0.5%. 

The structure of the paper is as follows. In Sec. II we develop the time-dependent mean-field scheme and derive 
coupled equations for the small-amplitude oscillations of the condensate and noncondensate components. Sec. Ill 
is devoted to spatially homogeneous systems. First we develop the perturbation scheme and hence we calculate to 
second order in the interaction coupling constant the equation of state of the system and the speed and damping 
rate of zero sound. In Sec. IV we apply the same perturbation scheme to harmonically trapped systems in the 
thermodynamic limit. We calculate the temperature dependence of the frequency shift of the low-energy collective 
modes and discuss the comparison with experiments. Finally, we show that at zero temperature our approach 
reproduces the hydrodynamic equations of superfluids. 



II. TIME-DEPENDENT MEAN-FIELD SCHEME 



Our starting point is the grand-canonical Hamiltonian of the system in the presence of an inhomogeneous external 
potential Vext{^)- In terms of the creation and annihilation particle field operators t) and tpir, t), the Hamiltonian 
takes the form 

+ dr^t(r,i)v,t(r,t)V'(r,i)V'(r,t) . (1) 

In the above equation we have assumed a point- like interaction between particles V{r — r') = gS{r — r'), with the 
coupling constant g given by the expression g — ATrh^a/m, to lowest order in the s-wave scattering length a. The 
equation of motion for the particle field operator follows directly from the Heisenberg equation and reads 
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The dynamic equations derived in this section correspond to the hnearized time-dependent Hartree-Fock-BogoUubov 
(TDHFB) approximation. This self-consistent mean-field scheme is based on the following prescriptions (we use the 
notations of Ref. 

a) V'(r,t) -*(r,t)+^(r,i) 
^r,t) = (^(r,t)) 
(^(r,i)) =0 

(^(r,i)^(r,i)) =TO(r,i) 

+m(r, t)iP^ (r, i)?/it (r, t) + m* (r, t)4>{r, t)iP{r, t) 

d) (^(r,0^(r,i)^(r,0> =0 
(^t(r,t)^(r,t)^(r,t)) =0 . 

The averages (...) in a), b) and d) are nonequilibrium averages, while time-independent equilibrium averages are 
indicated in this paper with the symbol (...)o- The prescription a) is the usual decomposition of the field operator 
into a condensate and a noncondensate component and defines the condensate wave function $(r, t). Prescription b) 
defines the normal, fi{v,t)^ and anomalous, TO(r, t), noncondensate particle densities. In terms of these quantities the 
interaction term in the Hamiltonian (^ quartic in the noncondensate components of '(/'(r, t) can be approximated using 
the factorization given by prescription c) . Finally, in prescription d) all averages of cubic products of noncondensate 
operators are set to zero. This is expected to be a good approximation for dilute systems. The inclusion of the triplet 
correlations {'>p'4>'4') and {'ip^'ipip) in a time-dependent self-consistent mean- field scheme is discussed in [^6[^. By using 
the above prescriptions one gets the following equation of motion for the condensate wave function 

in— $(r,t) - [ + Ve^tir) - m) $(r,i) +g|$(r,t)|2$(r,t) 

+ 2g$(r,t)n(r,i)+g$*(r,i)m(r,i) . (3) 

This equation includes the dynamic coupling between the condensate and the noncondensate particles. If we neglect 
these effects, n = rh = 0, equation (||) reduces to the usual Gross-Pitaevskii (GP) equation. 

We are interested in the small-amplitude oscillations of the condensate, which is only slightly displaced from its 
stationary value <l>o(r) = (V'(r))o 

$(r,i) = $o(r) +<5$(r,i) , (4) 

where (5$(r, t) is a small fluctuation. In the same way, we consider small fluctuations of the normal and anomalous 
particle densities 

n{r,t) = h°{r) +Sh{r,t) 

m{r,t) = m"{r) + 6rh{r,t) (5) 

around their equilibrium values ?T,°(r) — (■0'l'(r)V'(r))o and m°(r) = (-0(r)-0(r))o. 
The real wave function $o(r) satisfies the stationary equation |3l| ] 

+ Vext (r ) - /i + gno (r ) + 2gnO (r ) -f gmO (r) j $0 (r ) = , (6) 

where no(r) = |$o(r)P is the condensate density. The time-dependent equation for d^{r,t) is obtained by linearizing 
the equation of motion (||) 

z?i^<5$(r, t) = + K.t(r) - /i + 2gn(r)^ S^r, t) 

+ (gno(r) + gm°(r)) (5$*(r, t) + 2g$o(r)a(r, t) + g$o(r)<5m(r, t) , (7) 
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where we have introduced the total equiUbrium density n(r) = no(r) + n''(r). 

Both the stationary wave function $0 and the fluctuations 5$ depend through Eqs. and (0) on the normal and 
anomalous noncondensate particle densities, for which we need independent equations for their equilibrium values 
and time evolution. To this purpose it is convenient to express the noncondensate operators ifi, "ip^ in terms of 
quasiparticle operators ai, a| by means of the generalization to inhomogeneous systems of the Bogoliubov canonical 
transformations [P2| 



i 

Vit(r,t)=^(<(r)al(t)+«,(r)a,(t) 



(8) 



The normalization condition for the functions Ui (r) , Vi (r) , which ensures that the quasiparticle operators ai , aj satisfy 
Bose commutation relations, reads 



dr K(r)uj(r) - v*{r)vj{r)] = S,^ 



(9) 



The time evolution of n{r, t) and m(r, t) can be obtained from the Heisenberg equations for the products of quasipar- 
ticle operators a\{t)aj{t) and ai{t)aj{t) 



ih^^{alit)a,{t)) = {[alit)a,{t),H']) 
in^^{adt)aj{t)) = ([a,(t)a,(<),i/']> 



(10) 



In the above equations the commutators are calculated using the mean-field prescriptions a)-d) and the canonical 
transformation (H). The calculation can be easily done by noticing that in the Hamiltonian (|l|) only the terms 
quadratic and quartic in the noncondensate operators 'i/'j "0^ give non vanishing contributions, as we set to zero 
averages of single and cubic products of noncondensate operators. 

At equilibrium, we take the occupation of quasiparticle levels to be diagonal, {alaj)^ — Sijf^, while anomalous 

averages of quasiparticles are zero, {aiaj)o = 0. With these conditions, the stationary equations ( al{t)aj{t), H' )o = 
ai{t)aj{t), H' )o = yield the following coupled equations for the quasiparticle amplitudes Ui(r), Vi{r) 



Cuiiv) + [gno(r) -I- gfn° {v)]vi{v) = eiUi{r) , 
Cvt{r) + [gno(r) -I- gm°(r)]Mi(r) = -ejWi(r) 



where we have introduced the hermitian operator 



2rn 



+ Verf(r)-/i + 2gn(r) 



(11) 



(12) 



The coupled Eqs. ( pi] ) correspond to the static Hartree-Fock-Bogoliubov (HFB) equations as described in Ref. |jT^, and 
the Ci are the quasiparticle energies which fix the quasiparticle occupation numbers at equilibrium /? = [e'^*/'^^-^^ 1]^^- 
The equilibrium values of the normal and anomalous particle densities are written as 



(r)=5]{[|^..(r)|2^- |«.(r)n/°-f |z;,(r)|2} , 

i 

0(r) = E{2z.,(r)<(r)/° + z.,(r)<(r)} . 



(13) 



Out of equilibrium we define the following normal and anomalous quasiparticle distribution function 

/.,(<) = (a|(t)a,(t))-*y/0 , 



(14) 



in terms of which the fluctuations of n and m take the form 
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5n(r,t) = {K(r)",(r) + v* {r)v,{r)]f.,,{t) + u,{v)v,{v)g,,{t) + u*{v)v*{v)g*^{t)) , 
Srh{r,t) = J2 {2<(r)u,(r)/,,(t) + u,{r)uj{r)g,,{t) + <(r)t;*(r)g*.(t)} . 



(15) 



Up to linear terms in the fluctuations (5$, 6n and Srh, the equation of motion for fij gives the resuh [p3| 

X / dr $0 



(5$(i) ^UiW* + WiW* + + d^*{t) (^u^u* + v,v* + UiV*^ 



(16) 



2Sn{t) UiU* + ViV* j + S'm{t)viU* + Srh* {t)uiV* 



Analogously, for the time evolution of the anomalous quasiparticle distribution function gij , one obtains in the linear 
regime 



= (e, + e^9^J{t) + 2g(l + /f + /«) 



X dr 



Sm (^u*v* + v*u* + u*u*^ + 6'i>*{t) (^u*v* + v*v* + 



(17) 



+ g(l + /° + /°) / dv 25n{t)(^ 



Ui Vj + 



6fh{t)u*u* + Srh* {t)v* V* 



The first term on the r.h.s. of Eqs. (|T|) describes the free evolution of the quasiparticle states, corresponding 

to quasiparticle operators which evolve in time according to ai{t) = e~'^'^^*^^ai. The second and the third term 
describe, respectively, the coupling to the oscillations of the condensate and to the fluctuations Sh and 6rh of the 
normal and anomalous particle density. Above the Bose-Einstein transition temperature, where the system becomes 
normal and $ = m = 0, Eq. (16) corresponds to the linearized time-dependent Hartree-Fock (TDHF) equation. In the 
semiclassical limit TDHF is equivalent to the coUisionless Boltzmann equation for the particle distribution function 

il- 

In the framework of mean-field theories, coupled time-dependent equations for the condensate and noncondensate 
components of a Bose gas have been discussed by many authors. The TDHF scheme is discussed in |l^. Moreover, 
coupled equations of motion for the condensate and for correlation functions of pairs and triplets of noncondensate 
particles have been derived in ]T^ , and studied in the linear response regime in jiof . Similar kinetic equations 
were derived by Kane and Kadanoff ||35|] using the formalism of non-equilibrium Green's functions developed in ^6j . 
Recently, the approach of Kane and Kadanoff has been extended to deal with a trapped Bose-condensed gas in 1 37f . 

The stationary Eq. (^ for <f>o, with the normal and anomalous particle densities at equilibrium given by (|l^), and 
Eqs. ( |ll| ) for the quasiparticles correspond to the static self-consistent HFB approximation as reviewed by Griffin in 
1^ . These equations have been solved for harmonically trapped gases in |lj] . As it is well known |Q , in the case of 
homogeneous systems, the HFB excitation energies exhibit an unphysical gap at long wavelengths, which is fixed by 



the anomalous particle density: A = 2g^no\rh^\ (see Ref. |l2[). 

If one neglects the anomalous particle density, rh^ — 0, the HFB equations reduce to the so-called HFB-Popov 
approximation [p9lll2,Eq| , which is gapless and in the high temperature regime coincides with the Hartree-Fock scheme 



[ [40| . The HFB-Popov approximation has been applied by several authors to interacting bosons in harmonic traps, 
both to calculate the frequencies of the collective modes [|l^ and to study the thermodynamic properties of the system 
j^lj-^l .^ Gapless static mean-field approximations, alternative to HFB-Popov, have been put forward and discussed 

Finally, by neglecting both the normal and the anomalous particle density, nP = mP = 0, the HFB equations 
reduce to the Gross-Pitaevskii theory. From Eq. ^ one recovers the stationary GP equation, while Eqs. ( |TT1 ) follow 
from the linearization of the time-dependent GP equation. At very low temperatures, where effects arising from the 
depletion of the condensate are negligible, the Gross-Pitaevskii theory is well suited to describe dilute gases in traps. 
For these systems the linearized GP equation has been solved by many authors |^,^,|3| and successfully compared 
to experiments. 

The linearized TDHFB mean-field approximation is a closed set of self-consistent equations which describe the 
small oscillations of the system around the static HFB solution. The dynamic Eq. ^ describes the fiuctuation 5$ 
of the condensate around the stationary solution while Eqs. (|l5|)-(p7|) describe the small oscillations Sh, 6rh of 



5 



the normal and anomalous particle density around their equilibrium values nP, mP. Since the equations for the time 
evolution of (5$, dn and Sm are derived from the corresponding exact Heisenberg equations, it can be easily checked 
that the linearized TDHFB approach preserves important conservation laws, such as number of particles and energy 
conservation. This is a general feature of time-dependent mean-field approximations [^,^. Another important 
feature of linearized TDHFB is that, even though the quasiparticle energies obtained from Eqs. (|ll|) exhibit a gap at 
long wavelengths: ep ^ A as p ^ 0, the self-consistent solution of Eq. p) is gapless. In fact, one can show that this 



self-consistent solution satisfies the Hugenholtz-Pines theorem |47 



There are some questions one should address before embarking on the difficult task of a self-consistent solution of the 
linearized TDHFB equations. A first problem concerns the equilibrium anomalous density mP [see Eq. (H)], which 
is ultraviolet divergent. To second order in the interaction and for homogeneous systems the ultraviolet divergence 

is canceled by the renormalization of the coupling constant (see e.g. |48 ): g g (^1 + g"^ X]p "^/^'^) • How to 
include the renormalization of g in a self-consistent calculation and how to extend this renormalization procedure to 
inhomogeneous systems is still an open problem. Recently, Burnett and coworkers @,|l^ have shown that there is 
no need of renormalization of g if one uses, instead of a contact potential, an effective interaction, the many-body 
T-matrix, which includes self-consistently the effects arising from the anomalous density. Another problem concerns 
the gap exhibited by the quasiparticle energies in Eqs. (pT|). The self-consistent solution of these equations defines the 
equilibrium state of the system: it fixes the noncondensate densities and rhP through Eqs. ([T^), and the chemical 
potential /i and the condensate wavefunction <I>o through Eq. (^) . Even though the small oscillations around the static 
HFB solutions give rise to a gapless spectrum for the fluctuations of the condensate, properties such as the phonon 
velocity and their damping rate will be affected by an incorrect description of the system at equilibrium, originating 
from the unphysical gap A in the quasiparticle energies. In particular, if ksT ~ A, one expects a strong influence 
of the gap on the temperature dependence of these properties. In the present work we will not solve the linearized 
TDHFB equations self-consistently, instead, we solve them perturbativelly to order g^. By this method we avoid the 
problems mentioned above, in particular, the quasiparticle states in (11) are properly described by Bogoliubov theory, 
which is gapless. 

Another point which deserves some comments is the Kohn mode (dipole mode). As it is well known, in the presence 
of harmonic confinement the center of mass degrees of freedom separate from all other degrees of freedom, giving rise 
to a collective mode of the system, the Kohn mode, in which the equilibrium density profile oscillates rigidly at the 
trap frequency. The linearized TDHFB equations obtained in this section do not describe this mode, as they do not 
account for the motion of the center of mass. In fact, these equations correctly describe excitations for which the 
center of mass is at rest, and we will restrict our analysis to this class of excitations. The Kohn mode is associated 
with broken translational symmetry and is often referred to as a "spurious" mode. For a discussion of spurious states 
and their appearance in linearized time-dependent mean- field theories see Ref. p6|. 



III. SPATIALLY HOMOGENEOUS SYSTEM 



In this section we will develop, starting from the dynamic equations for the condensate and noncondensate compo- 
nents derived in the previous section, a perturbation scheme for the elementary excitations in a homogenous system. 
Explicit results for the temperature dependence of the chemical potential, damping rate and speed of zero sound are 
given in the limit a^no <C 1, where uq = ho(T') is the equilibrium value of the condensate density at temperature T. 
At zero temperature, the calculation presented here is equivalent to the Beliaev second-order approximation of the 
single-particle Green's function p9[. In the high temperature regime, ksT 3> /i, our approach corresponds to the 



finite-temperature extension of the Beliaev approximation recently employed in Refs. |28 2^ (the former reference 
gives also a systematic account of the various perturbation schemes for a uniform Bosc gas). The perturbation ex- 
pansion carried out in the present work holds to order {a^noY^'^ for any temperature in the condensed phase, except 
clearly very close to the transition, where the time-dependent mean-field equations we used as starting point break 
down. 

In a spatially homogeneous system the condensate wave function at equilibrium is constant, $o(r) = y/nQ- The 
stationary equation (^) reads then 

^l = gno + 2gn° + gm° , (18) 

and represents the equation of state of the system, which fixes the chemical potential as a function of the condensate 
density no and the temperature T. By using the above result for the chemical potential, equation (|^) for the 
fluctuations of the condensate becomes 



6 



d ( ?i^V^ \ 

5$(r, t) = + gno - gm" (5$(r, i) + (gno + gm°) <5$* (r, i) 

of \ 2m / 

+ 2gVr^5n(r, i) + g^A^(5m(r, i) . (19) 

In the above equation the terms involving the equihbrium anomalous density m° account for the coupling between 
the fluctuations of the condensate and the static distribution of noncondensate atoms, while the terms containing 8h 
and (5m describe the dynamic coupling between the condensate and the fluctuations of the noncondensate component. 



A. Perturbation scheme 



The perturbation scheme applied to Eqs. (|18D, (|19| ) consists in treating the terms which give the static and dynamic 
coupling to the noncondensate component to second order in g. It means that the static and fluctuating parts of the 
normal and anomalous density need to be calculated only to first order in g. To accomplish this task one must retain 
in the equations for the quasiparticles (|ll|), ( p^ ) and ( p7| ) only the terms which describe the coupling to the condensate 
and neglect all terms arising from the coupling to the noncondensate component. 

Let us suppose that the condensate oscillates with frequency w and wave vector q/?i 



Furthermore, the quasiparticle amplitudes be described by plane-wave functions 



V V" 



-p(r) = 



(20) 



(21) 



To first order in g the quasiparticle equations (11) become 



(p /2m + gna)up 
(p^/2m + gno)wp - 



gnoVp 



(22) 



These coupled equations coincide with the well-known Bogoliubov equations for the real quasiparticle amplitudes Up, 



which satisfy the following relations 



Up = 1 + f p 



i4 + g'niy/' + ' 
2e„ 



Up Dp 



with the quasiparticle energy Cp given by the Bogoliubov spectrum 



{(4 



g"0 



1/2 



(23) 



(24) 



being = /2m the free-particle energy. Notice that, by employing the equation of state ([l^), the full HFB equations 
( pi] ) would coincide with the matrix equation ( |2^ ) apart from a term gmg. This term would appear with the minus 
sign in the diagonal term and with the plus sign in the off-diagonal term, and is responsible for the gap in ep as p — s- 0. 
Since we use the Bogoliubov result (p^, we avoid the problem of the gap in the quasiparticle spectrum. 

In the same approximation one must neglect in Eqs. ( p^ ) and (17) the terms which describe the coupling to the 
fluctuations 5n and 5rh of the noncondensate component. Due to the coupling to the condensate, which acts as a 
time-dependent external drive, the only elements of the matrices /p,p', ffp.p' and 5* p, which oscillate at the frequency 
Lo are given by 



/p,q+p(^) — 



/p f\ 



|q+p| 



(l5$l + 6^2) (2MpU|q+p| + 2i;pW|q+p| + WpU|q+p| + Mp^^|q+p|] 
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5p,q-p(w) = g 



y/n^ 1 + /p + /|°q_p| 



((5$l - 6^2) {UpU\q-p\ - VpV\ 



y/V huj - (ep + e|q_p|) + iO 

+ (5$2) (2MpW|q_p| + 2WpU|q_p| + UpM|q_p| + fpW|q-p|) 



q-p 



(25) 



( ] = ^ J- + /" + /|q+p| 

5p,-q-pM - g ^ ^ ^ . 



((5<I>1 - ^$2) (MpM|q+p| - fpt^|q+p|] 



- ((5i>l + (5$2) (2Mpi;|q+p| + 2t;p-U|q+p| + UpU|q+pl + VpV\c^+p\ 

In the above equations the frequency w has been chosen with an infinitesimally small component on the positive 
imaginary axis. 

As it is well known (see e.g. Ref. to treat consistently to order g^ the properties of a Bose-condensed gas one 

must include the proper renormalization of the coupling constant g — > g ^1 + g-^- m/p^^ . This renormalization 

of g is crucial because it cancels exactly the large-p ultraviolet divergence exhibited by the equilibrium anomalous 
density ?7i°. In fact, by using the renormalized value of g, the term gno + gm° present in Eqs. (|l8|), (foh becomes 



gno + gm -> gno + g nQ— I - 

p 



m 1 + 2/p"^ 



2e. 



(26) 



and is well behaved at large p. 

To order g^, Eq. (|l^) for the fluctuations of the condensate can be finally written in the form 



hu) ((5$i + 6^2) = TT— (^*1 - <^*2) 

2 m 



Eii(q,w) - Eii(q,-w) 



((5$i + (5$2) 



Sii(q,cj) + Eii(q,-[j) 



^ + 2gno j ((5$i + ^$2 
I]ii(q,cj) + Sii(q, -w) 



Si2(q,t^) 



{d<^i - (5$2) , 



{6^1 - 5^2) 



(27) 



+ I]i2(q,t^) 



((5<I>i + (5$2) 



In the above equations the self-energies I]ii(q, w) and Si2(q, w) are proportional to g^. They are obtained from Eq. 
(19) through the expressions (15), which give 6h and dm in terms of the matrices /p,p' and gp,p', and through Eqs. 
(25) and (|2q). After some straightforward algebra one finds for the self-energy Sn 



2gpAfc + SCpAk + 4Bpgfc + 4CpCk 2ApBk + SCpg^ + iApA^ + 4CpCfc 
-f (ep — gfe) -I- iO — (ep — gfe) -f- iO 



no 



(28) 



p p 

2ApAk + SCpAk + 4ApBk + iCpCk 2BpBk + 8CpBk + ABpA^ + 4CpCfc 



huj — (ep -I- efc) -I- iO 



?iw -I- (ep + ek) + »0 



where we have introduced the notations: k = q + p, Ap = Up, Bp = Vp and Cp — UpVp. Analogously for S12 one has 



I^12(q, .) = g^no^ E f S - ^) + E - 

P \ / P 



4CpBfe -f 4Cp^fe + ABpBu + eCpCfc 4CpAfc + ACpBu + AApAu + (SCpCk 



huj + (ep — efc) -I- zO 



?ia; — (cp — Efc) + iO 
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g'-o^^(l + /p+/") (29) 
p 

/ iCpBk + iCpAk + AApBk + GCpCk ^CpB^ + ACpA^ + iBpAk + 6CpCk 



\ huj - {tp + Cfc) + zO hoj + {tp + efc) + iO 

The above expressions for Sn and S12 coincide with the second-order self-energies exphcitly calculated at finite 
temperature by Shi and Griffin using diagrammatic methods . At zero temperature they correspond to Beliaev's 
results [Q, while in the high-temperature regime, fc^T ^ gno, they have been recently discussed by Fedichev and 
Shlyapnikov 

By neglecting in ( p7| ) the terms proportional to the self-energies, one is left with the equations for the fluctuations 
of the condensate to first order in g. These equations coincide with the quasiparticle equations (p2|). The solution is 
then given by 5$i(q) = Uq and J$2(q) = Vq, with Uq and Vq given by (p3|). The excitation energy is given by the 
Bogoliubov spectrum (|24|). To order g^ one writes 

huj = Eq + SE-ij , (30) 

where SE is the real part of the frequency shift and 7 is the damping rate. It is straightforward to obtain the 
second-order correction to huj from Eqs. (p7|), one finds 

SE~ij = i:ii{ci,eq)ul + 2J:i2{ci,eq)uqVq + Eii{q,-eq)v^ , (31) 

where the self-energies have been calculated for huj — Cq. After some algebra one can cast Eq. ( ^l[ ) in the more 
convenient form 

5E-i-i= {uq + Vqf g^no^ ^ " " ^9)^ 8"*° 

p ^ 

+ 2g'lj;(l + />/,») ( , °X'\-.n - ^/f'l^n ) ' 

y \eq - [ip + ek) + Cq + (cp + Cfe) + lOy 

where we use k = q -f p and we have introduced the matrices 



ApM = — ^ [{Uq + Vq){2UpUk + 2VpVk + VpUk + UpVk) + [Uq - Vq){VpUk - UpVk)] 



Bp,k = — ^ [{Uq + Vq){2UpVk + 2VpUk + UpUk + VpVk) + {Uq - Vq){UpUk - VpVk)] , (33) 



Bp,k = [{Uq + Vq){2UpVk + 2VpUk + UpUk + VpVk) - {Uq ~ Vq){UpUk - VpVk)] ■ 

Result ( |3^ ) gives the correction to the Bogoliubov elementary excitation energy tq. The first and second term 
on the r.h.s. of (^2|) arise, respectively, from the renormalization of g and from the equilibrium anomalous density 
•fhP . The last two terms arise, instead, from the dynamic coupling between the condensate and the noncondensate 
component. The real part of the r.h.s. of (^) gives the frequency shift 5E^ while the imaginary part gives the damping 
rate 7. Notice that, concerning the damping rate, result (|3^ ) coincides with the calculation carried out within the 
Popov approximation [Eq. (39) of where the condition fhP = was assumed. However, as discussed in [ p^ , 

the frequency shift 5E is not given correctly to order g^ by the Popov approximation. In fact, the static anomalous 
density fhP and the renormalized coupling constant contribute to the real part of Eq. (32). 



B. Equation of state 

Let us first discuss the equation of state (p^). By calculating vP and fhP using the equilibrium expressions (|l 
with the first-order quasiparticle amplitudes and energies given by (p3|) and (|2^), one finds to order g^ 



= gno + 2gn!^-Kgno(a3no)i/2^(T) . (34) 
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In the above equation n!^ 



C(3/2)Ay^ ~ 2.612Ajn'^ is the noncondensate density of an ideal gas, which is fixed by 



the thermal wavelength At = ^iy^Tr/mfceT. Moreover, H{t) is a dimensionless function of the reduced temperature 
T = ksT/gno given by 



H{r) = 



40 



32 



dx- 



- 1 



-2u- 1 



1- 



(35) 



where we have introduced the quantity u = Vl + r^x^. Result (|34| ) gives the chemical potential as a function of the 
equilibrium condensate density uq and the temperature T to second order in g. It coincides with the result of Shi 
and Griffin [Eq. (7.9) of |2^]. Notice that the sum of the first two terms on the r.h.s. of (jsj) corresponds to the 
chemical potential calculated to first order in g. In Fig. 1 the dimensionless function H{t) is plotted as a function of 
the reduced temperature r. 

At low temperatures, r ^ 1, the function H{t) can be expanded as 



Hir) 



40 



32C(3/2)t 



3/2 



27r3/2 
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In the same regime of temperatures, the condensate density is given in terms of the total density n = hq 
following relation 



(36) 



h by the 



no 



27r3/2 



30F 



30 



(37) 



valid to order g^. In the above expression the first term in brackets corresponds to the quantum depletion of the 
condensate, while the other two terms account for the thermal depletion caused by phonon-type excitations. By using 
Eqs. ( ^ ) and (|3^), one gets the following result for the low-temperature behavior of the chemical potential in terms 
of the density n 



l + (a%)^/' 



/ 32 2^7/2 



= 0) + 



V30F 
60 nofi 



15 



3^3 



(38) 



where /x(T = 0) = gn[l + 32(a'^no)^^^/3\/7r] is the value of the chemical potential at zero temperature [Q, and 
cb = yJgriQ/m is the Bogoliubov velocity of sound. The term in ( ^8|) coincides with the result obtained from the 
thermodynamic relation /i = (dF / dN)v.T, where F is the low-temperature expansion of the free energy of a Bose gas 



At high temperatures, r 3> 1, the function H{t) yields the asymptotic result H{t) 
potential in this regime of temperatures one gets 

/.t ~ gno + 2gn§. - l2V^{a^mf'^kBT , 

which coincides with the result obtained by Popov |39,Ba,E^. 



-12-y7rT, and for the chemical 
(39) 



C. Zero sound: damping rate and velocity of propagation 

Zero sound, or Bogoliubov sound, is a collective oscillation of the system in the coUisionless regime, for which the 
restoring force acting on a given particle comes from the mean-field created by the other particles. At zero temperature, 
zero sound coincides with ordinary sound and the velocity of propagation c is fixed by the compressibility of the system 
mc^ = {dP/dn)T=o- At finite temperature one can have both zero sound and hydrodynamic modes, depending on 
the temperature and the wavelength of the excitation. In this case the velocity of zero sound can not be fixed only by 
thermodynamic relations. For an exhaustive discussion of coUisionless and hydrodynamic collective modes we refer 
to the books [|jl^]. 

To first order in g, the excitation energy of the zero-sound mode is given by the long-wavelength limit of the 
Bogoliubov dispersion relation (p^, which corresponds to phonons eg = cbQ propagating with the velocity cb = 
a/ gno / ni fixed by the condensate density. Starting from Eq. (^^, we will explicitly calculate the damping rate and 
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the temperature dependence of the speed of zero sound to order g^. The damping of phonons in homogeneous systems 
has been recently investigated by Liu using functional-integration methods, and by Pitaevskii and Stringari |2^] 
by means of perturbation theory. The temperature dependence of the speed of zero sound has been investigated by 
Payne and Griffin within the framework of the dielectric formalism, and by Shi and Griffin |2|] and Fedichcv and 
Shlyapnikov using diagrammatic methods. 



1. Damping of zero sauna 



The calculation of the damping of phonons from Eq. (|3^) has been already carried out in Here we will only 
review the main results. 

In the quantum regime, ^ ksT, the damping is governed by the Beliaev mechanism, in which a phonon decays 
into a pair of excitations. This mechanism is described in Eq. (|32| ) by the imaginary part of the term containing the 
matrices B and B. The damping rate is given by 

^ = ¥^ . (40) 

This result has been first obtained by Beliaev using diagrammatic techniques |Q . 

In the thermal regime, Cg ksT, the phonon decay is dominated by a different damping process, in which a phonon 
with energy Cq is absorbed by a thermal excitation with energy ep and is turned into another thermal excitation with 
energy e|q+p|. This mechanism is known as Landau damping and is described in Eq. ( |3^ ) by the imaginary part of 
the term involving the matrix A. The result is given by (see Refs. [ pO|jl^ ) 

^ = [a^nof/^F{T) , (41) 
esq 

In the above equation r = ksT /gn^ is the reduced temperature and Fir) is the dimensionless function 



(42) 



where u is defined as in (|35|). 

For temperatures r <C 1 the function F takes the asymptotic limit F — > 37r^/^r''/5 and one recovers the known 
result for the phonon damping j5^-|5^,H 



7 _ 37r3 {keTY 



cbQ 40 mnoh c 



3^5 



(43) 



In the opposite regime of temperatures, r 3> 1, one finds the asymptotic limit F 37r^/'^r/4, and the damping 
coefficient is given by 

7 SnksTa ^^^^ 



CbQ 8 he 



B 



The damping of phonons in this regime of temperatures has been first investigated by Szepfalusy and Kondor | |5q |. 

In Fig. 2 the dimensionless function F{t) is plotted as a function of t together with its asymptotic behaviour both 
at small and large t's. One can see that F departs rather soon from the low-temperature t** behaviour, while it 
approaches the high-temperature linear law very slowly. 



2. Velocity of zero sound 

Differently from the calculation of the damping rates, all terms on the r.h.s. of (^2|) contribute to the frequency 
shift SE. The first two terms, which involve mo and the renormalization of the coupling constant, are referred to as 
static terms. Concerning the other two terms, for excitations with <?C gno, the relevant region in the calculation of 
the real part of the term which containes the matrices B and B corresponds to energies Cq. This term is referred 
to as non-resonant term. On the contrary, the resonance region gives an important contribution to the real part of 
the term involving the matrix A. This term is referred to as resonant term. 
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Let us start by analysing the non-resonant term. The contribution of this term to the energy shift 5E can be 
written as 



2e„e 



2e„ 



-^{2upVk + 2vpUk + UpUk + VpVkf 



Cp + Cfe 



,0 , , N , (UpV^-VpVkf 

[2upVk + 2vpUk + UpUk + VpVk) + — — -, — — — — 



(45) 



where k = q + p. The above result is valid for any excitation energy Cg, and is not limited to the long- wavelength 
regime €q <^ gf^o■ In the phonon regime, result (|45|) can be simplified and one gets 



esq 



p 



0\2 



463 



1 2 1 V- 



1^1 + 2/0 



2e„ 



2epek 



(46) 



The contribution to 6E from the resonant term can be calculated in the same way and one gets the general result 



f 



Ep - Cfe Jp Jk ^p) 



2e„efe 



1 fO _ fO 

, 2 ^ Jp Jfc 

V ^ eq + ep - ek 



[2upUk -f 2?;pWfc -I- VpUk + UpVk ) - 2— — 



{2upUk + 2vpVk + VpUk + UpVk) 



(ep - efe)2 



(47) 



In the limit <C gJio the above expression reduces to 



esq ~ 2V 



^ dep 



^g'"o-^E/p- 



2tpek 



1 - 



CB 



cb — cos 9 dep /dp 



2 gnpeO 

3 e2 



(48) 



where is the angle the momentum p forms with the direction of q. 

Notice that the last terms on the r.h.s. of ( p6| ) and ji^ ) are equal and opposite in sign, thus they cancel in the 
sum 5Ept + 6Enr- Moreover, the second term on the r.h.s. of ([iq), which is ultraviolet divergent, is canceled by a 
corresponding term arising from the contribution to SE of the static terms. In conclusion, the final result for the 
energy shift SE in the phonon regime is well behaved and proportional to e^. One finds 



CBq 2V 



—y 

2V ^ 



dfl 




dep 


\ep 


1 m 


1 




" 2^ 



de" 



CB 



cb — COS 6dep/ dp 



inoe" 



2 gnpe" 
3^ 



3 £3 



3 el 



(49) 



The r.h.s. of the above equation gives the correction to the Bogoliubov velocity of zero sound. By rearranging the 
integrals over momenta, one gets the relevant result 



CB 



1 + ia^no)'^^G{T) 



(50) 
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Here G(r) is the following dimensionless function of the reduced temperature t = fcsT/gno 
28 



G(r) 



30F 



32 f°° ^ 1 Vti^ 5 
— r / ax— 

Jo e' 



3u 



1 



3u\/u + 1 



^ 3 (2m+ l)^(u- 1) 



u 6(u + 1) 
1 



dx- 



1 



Vm + 1 



log 



V2m - VuTT ' 



\/2m + + 1 



(51) 



where u is defined as in (|35|). 

It is interesting to study Eq. ( pO[ ) in particular regimes of temperature. At zero temperature, the function G takes the 
value: G{t = 0) = 28/(3-^/7?), and no is related to the total density n by the expression: no = n[l — 8/(3-y7r)(a^no)"^^^], 
which accounts for the quantum depletion. The result for the sound velocity is: 



e(r = 0) = 

V m 



1 + ^(a'^no 



,1/2 



(52) 



The above result, which was first found by Beliaev coincides with the one obtained from the thermodynamic 
relation c(T = 0) = [n(5/i(r = 0)/dn)/m]^/'^, where fj.{T = 0) is given in (|8|). 
At low temperatures, t ^ 1, one finds the following expansion of the G function 



G(r) 



28 



r3/2 



37rV2 



r^log(l/r2) 



(53) 



In this regime of temperatures the condensate density is given by the expression (|3j) in terms of the density n, and 
the velocity of zero sound turns out to be 



<^ - 0) + ^ JMJl iog[^2^V(fc^r)2] 



(54) 



This result was first obtained by Andreev and Khalatnikov using kinetic equations, and later by Ma et al. |57| ] 
within the framework of the dielectric formalism. 

Finally, in the high temperature regime r ^ 1, the function G is linear in r: G(r) — > G(oo)r, with the numerical 
coeflticent G{oo) given by the following expression 



G(oo) ^(9%/2-28) + ^ / dx- 
3 Jo 



3a;2 - 4 , / \/2 - Jx{l + x) 
log ' 



Vl"a;2(l + a;) \^+^x{l + x) 
For the speed of zero sound in this regime of temperatures one gets the result 

G(oo) ksTa 



-7.4 . 



CB 



20F h 
2.1 agrees with the finding of 



(55) 



(56) 



while is about a factor 6 larger 



where the numerical coefficient G(oo)/(2-y/7r) 
than the one calculated in ||28|| . 

The proper description of the cross-over between the low and high-temperature regime is provided by Eq. ( |5C 
The dimensionless function G(t) is plotted in Fig. 3. In the experiments on trapped gases the gas parameter in 
the center of the trap is typically a'^no ^ 10~^-10~''. For temperatures of the order of the chemical potential, which 
means r ~ 1, the correction to the Bogoliubov speed of sound amounts to about 2-5%. 



IV. SPATIALLY INHOMOGENEOUS SYSTEM 



In this section we generalize the perturbation scheme developed for a homogeneous Bose-condensed gas to the case 
of inhomogenous systems trapped by a harmonic confining potential 



(57) 



The relevant length scale associated to the external potential (B^) is the harmonic oscillator length defined as 
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where ujho = {t^xi^yi^zY^^ is the geometric average of the oscillator frequencies. The length scale aho gives the average 
width of the Gaussian which describes the ground state of non-interacting particles in the harmonic potential (pj). 
The shape of the potential V^xt fixes the symmetry of the problem. So far all experiments on trapped Bose gases 
have been realized using axially symmetric traps. In this case there are only two distinct oscillator frequencies: 
u;± = Wx = ojy and uJz- The ratio between the axial and radial frequencies, A — uJz/lo^, fixes the asymmetry of the 
trap. For A < 1 the trap is cigar shaped, whereas for A > 1 is disk shaped and A = 1 refers to spherically symmetric 
traps. 

In our analysis we only consider systems with repulsive interactions (a > 0) in the thermodynamic limit. As 
extensively discussed in for harmonically trapped Bose systems the thermodynamic limit is obtained by letting 
the total number of trapped particles N ^ oo and cuho 0, while keeping the product Nujf^^ constant. With this 
definition the Bose-Einstein transition temperature fc^T^ = ?iti;/jo(A^/C(3))^/^ is well defined in the thermodynamic 
limit. 

In the thermodynamic limit, the condition Ni^{T)a/ ^ 1, which ensures the validity of the Thomas-Fermi (TF) 
approximation for the condensate with occupation number A'o, is always guaranteed below the transition temperature. 
In the TF approximation one neglects the quantum-pressure term proportional to V^$o(r) in the stationary equation 
and the equilibrium profile of the condensate density is fixed by the following equation 

gno(r)=/^-Fext(r)-2gnO(r)-gmO(r) , (59) 

in the central region of the trap where the r.h.s. of the above equation is positive, whereas outside this region one has 



no(r) = 0. The chemical potential in Eq. (59) is fixed by the normalization condition J dr no(r) = iVo(T), with Nq{T) 
the equilibrium condensate occupation number at temperature T. To lowest order in the interaction, the profile of 
the condensate density has the form of the inverted parabola 

nTF{r)^g-'biTF{No)~Vext{r)] , (60) 

where 

hujiio f IBNqq^ 



/XTF(iVo) = ^ ( ] (61) 



is the temperature dependent TF chemical potential. Moreover, in the thermodynamic limit, one can show [p2P] 
that the equilibrium properties of the system can be expressed in terms of two parameters: the reduced temperature 
t = T jT^ and the interaction parameter ry defined as the ratio 

Mtf(^) 

between the TF chemical potential at T = and the transition temperature. 

The time-dependent equation (Q) for the fiuctuations of the condensate only needs to be solved in the region where 
no(r) 7^ 0, according to Eq. ([59|). One finds 



ih—5<^{Y,t) 



gno(r) - gm"(r)^ ^<I>(r, t) + (gno(r) + gm°(r)) (5$*(r, t) 



dt " ' \ 2m 

+ 2g$o(r)(5n(r,i)+g$o(r)(5m(r,i) . (63) 



For a trapped system in the thermodynamic limit the above equations (59) and ( p3| ) replace respectively Eqs. ( p^ ) 
and (pj|), holding for a homogeneous system. 

We are interested in the lowest-lying collective modes of the system with excitation energy fno ^ /i. To lowest 
order in the interaction these modes are the solution of the following coupled equations 

z;i^[<5$(r,i) + <5$*(r,<)] ^ [5$(r, i) - ^<i>*(r, i)] , 

in^[<5$(r,i)-<5$*(r,t)] =2gnTF(r)[(5$(r,t)+5$*(r,i)] . (64) 
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These equations are obtained from ( |63D by neglecting all coupling terms to noncondensate particles and neglecting 
also the term proportional to V^((5$ + (5$*), which is of higher order for the low-lying modes we are considering [^ . 

The oscillating solution defined by i5<i>(r,t) = S^i{r)e^^'^* , S^*(r,t) — (5$^(r)e^*'^*, with the Fourier components 
fixed by the relations 



(SitlT) - i*»(t)) = y^S^g^ Xo(r) , (65) 

reduces the coupled equations (^4|) to the following equation for the function xo(r) 

mio\o (r) + V [gUTf (r) Vxo (r)] = . (66) 

The normalization condition J dr (|(5$5P ^ l<^*J'2p) = 1 satisfied by the Fourier components and (5$2i implies the 
normalization condition J dr XoXo 1 on the function Xo(r). 

Equation ( |66| ) was first derived at T = by Stringari using the hydrodynamic theory of superfluids, and it has 
then been studied by many authors [^,^. For spherically symmetric traps the excitation energies hui = exF obey 
the dispersion law ||6| 

1 /2 

erpinrj) — fiuJho {"^nf. + 2nrl + 3nr + I) , (67) 

where and / are respectively the radial and the angular momentum quantum numbers. In the case of axially 
symmetric traps, analytic results for the excitation energies are obtained for the m = low and high mode B 



CTFim = 0)lm = hoj^ (^2 + ^X^ T ^V9A4-16A2 + 16^ 



1/2 

(68) 



and for the surface modes of the form Xm T-I^'le'""^, for which one has 



eTF{m) — hu!±y^\m\ . (69) 

A general feature of Eq. (^6|), which is explicitly reflected in the above results for exF, is that the excitation energies 
do not depend on interaction and are proportional to the oscillator frequencies of the harmonic potential. At finite 
temperature, where griTF — ^J'TF[No{T)] — Vext, this fact implies that exF does not depend on temperature either. 
This is an important difference with respect to the homogeneous case where the corresponding excitations have the 
dispersion law eg — ^ griQ/m q, and depend on temperature through the condensate density. The behavior exhibited 
in the harmonic trap is well understood if one notes that the values of q are fixed by the boundary and vary as 1/R, 
where R is the size of the condensate. In the Thomas-Fermi limit, R ~ \/ fir f / "fni-^ho and the radius R explicitly 
depends on the ch emical p otential. On the other hand, the sound velocity is also fixed by the value of the chemical 
potential: ~ a/ ij,tf/iti. One finally finds that in the product esq the chemical potential cancels out, so that the 
collective frequency is proportional to the oscillator frequency ujho- 

Provided that finite size effects can be neglected, the explicit dependence of the collective frequency on the inter- 
action parameter r] as well as on the reduced temperature t = T/T^ arises due to quantum and thermal fluctuations 
which are of order g^ . These fluctuations have the same physical origin as the corrections to the Bogoliubov speed of 
sound given in the homogeneous case by result (|50|). The difference, however, is that in the case of harmonic traps, 
beyond mean-field effects give corrections to a collective frequency which is fixed only by the oscillator frequency: a 
much better situation from the experimental point of view. In the following part of this section we will explicitly 
calculate the effects of quantum and thermal fluctuations on the frequencies of the lowest compressional and surface 
modes. 



A. Perturbation scheme 



The perturbation scheme we employ for trapped systems follows the same lines as the one developed in the homo- 
geneous case. However, there are two important differences: first of all the quasiparticle states are not exact plane 
waves, secondly the condensate density is not fixed by a single parameter but depends on position. 
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Concerning thequasiparticle states, we make use of the local density (semiclassical) approximation which amounts 



to setting lira 



e^-^-W ^ = ^e^'^'W ^ (70) 



where is a large volume contaning the system and the real functions Ui, Vi satisfy the normalization condition 
(r) — vf{v) — 1. The factor e*'''''^''^ represents the rapidly varying part of the functions Ui and w^, while the functions 
Mi, Vi are assumed to be smooth functions of the position. The phase ^i{v), which is also assumed to be a smooth 
function of r, characterizes the local impulse p = hVipi of the quasiparticle. When summations over quasiparticle 
states are involved, these are replaced in the semiclassical approximation by sums over momenta, ... — > 
and Ui{v) — > Up(r), ^^(r) — s- Vp{v), where the functions Up{v) and Wp(r) are given, in the region of the condensate, by 
the following expressions 

2, , 2, ,_{el{v)+g^nl^{v)Y'^+ep{v) 



2ep(r) 



'^.W-pW - , (71) 



2ep(r) 

where the position-dependent quasiparticle energies ep(r) are given by 



ep{v) = ([el + gnTF{v)Y-[gnTF{v)\')"^ . (72) 



For each position r the above equations coincide with the Bogoliubov expressions (Pg), ( p4| ) with a local condensate 
density given by the TF density profile titf (r) defined in (|6^) . The semiclassical approximation for the excited states 
of a trapped Bose gas has been extensively used in the theoretical study of the thermodynamic properties of the 
system 1 4l|, p2| , |6l| . It gives a very good description of the system for temperatures ksT hojho, but is also vahd at 
T = if the relevant energies in the summation over excited states are much larger than the oscillator energy huho 
p3[ . For large systems the oscillator energy is the smallest energy scale and vanishes in the thermodynamic limit, as 
a consequence, in this limit, the semiclassical approximation becomes a rigorous treatment. 

The equilibrium noncondensate densities tiP(y) and m''(r) are readily calculated employing the semiclassical ap- 
proximation. One obtains 



p 



= -^E"'W^'W(l + 2/°) = ^E"p(^HW [l + 2/p"(r)] , (73) 



I p 



where /p (r) = (e'^p^"")/'^^'^ — 1)"-'^ is the local equilibrium quasiparticle distribution function. 

By inserting in Eq. (^9|) the above expressions for the noncondensate densities and using the renormalization (|2^ ) 
of the coupling constant, one obtains the following result for the profile of the condensate density valid to order 

gno(r) = gnTj=^(r) +3^- gnTF{r)[a^nTF{r)]^^^ H[T{r)] , (74) 

where H{t) is the dimensionlcss function ( |35| ) of the local reduced temperature T(r) = ksT/griTFi^)- In the above 
equation Sfi = fj,— iiTFiNo) — 2gn^, with — (^{S/2)X^^ as in (^), is the shift in the chemical potential corresponding 
to the change in the condensate density profile. 

The application of the semiclassical approximation to the last two terms on the r.h.s. of Eq. (|6^), which describe 
the dynamic coupling to the noncondensate particles, needs a careful treatment. ThuSjfor the moment, we calculate 
them in terms of the Ui and Vi functions. Similarly to the homogeneous case [see Eqs. (|2^)] , one must neglect in Eqs. 
(p^), ( p^ ) the terms proportional to Sh and Srh. In Fourier space, one finds for the components of the matrices /y, 
gij and g*j oscillating at the frequency ui 

M.) = g^J^f^fy^ / dr ^0 [iS^. - S<!>,)i..u* - ..«;) 
+ ((5$i + S^2)i2uiU* + 2v^v* + v^u* + UiV*)] , 
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5y (t^) 



9h (^) 



J dv'fo[{6^i-S<P2)iu*u*-v*v*) 
*v* +2v*u*+u*u*+v*v*)] , 



hu> — {ti + ej) 
+ + 54>2)(2<w* + 2w>* + it>* + v*v*)\ 
^ 1 + /° + /° 



(75) 



where we have taken ^$(r, i) = (5$i(r) e~*'^* and (5<i>*(r,t) = ^$2(r)e~*"*. In Eqs. ( [75| ) we have neglected in the 
expressions for gij{uj) and g*j{i-o) the small imaginary part of the frequency co. As discussed in Sec. III-C, this 
imaginary contribution is responsible for the Beliaev damping. However, due to discretization of levels, this damping 
mechanism is not effective for the low-lying modes we are investigating. 

To order g^ the equation for the low-lying oscillations of the condensate can be written in the form 



huj{d^i + 6^2) = - 



2m 



■(^$1 - (5$2) 



?iw((5<i>i 



- 2gTO°((5$i - 5*2) + 2gV^^(v>j - u*Vj)f,j(uj) 
(5$2) = HgriTF + S^i){S^i + 6^2) 



(76) 



2gnTF 



((5$i + 6^2) 



p 

+ 2gVHi^^(2u>j -I- 2v*Vj + v*Uj + u*vj)f,j{u;) 
+ gy/riTF |^(2uiWj -I- 2viUj + UiUj + ViVj)gij{u)) 

ij 

+ (2<w* + 2v*u* + u*u* + v*v*)gl^{oj) 

In the above equations one can recognize the terms arising from the dynamic coupling between the condensate and 
the noncondensate component, which contain fij(uj), gij{i-u) and g*j{uj), the terms arising from the coupling to the 
static anomalous density rh'^ and from the rcnormalization of g, and, finally, the terms proportional to and H{t) 
which come from the change in the density profile of the condensate. The last terms have no counterpart in the 
homogeneous case. In Eqs. (|7|) we have neglected, as in Eqs. (|64[), the term proportional to W^{6^i + 6^2)- 

Following the analysis carried out in the homogeneous case, we write the excitation energy as hu = eTP + 5E — 17. 
From Eqs. (JTq), by treating the corrections to Eqs. (64) as small perturbations, one gets the result 



6E~ij^ JdrgriTF g^"^^ - {a^nTFy^'^H[T{r)] 

- /drgmV$;-^<i>^P + 4g^E(/°-/i)— ^^0 
J ^ ^ ETF + (e« - ej ) +»0 

+ 2g^5:(l + /f + /j') ' 



(77) 



eTF-{£i + <^j) (TF + {^i + ^j) 
"J \ / 

holding for the low-lying modes with e^F Notice that the shift dpi of the chemical potential does not enter result 

([77I). In fact, in the Thomas-Fermi limit, the excitation frequencies obtained from Eq. (|6^ ) do not depend on the 
value of /i. In Eq. ([77|), the matrix elements Aij, Bij and _By are defined, in analogy to the homogeneous case, as 

1 
2 
1 
2 



Aij j dr ((5$? -f <5$^)(2m,u* -I- 2v,v* + v,u* + mv*) -t- ((5$? - 5<^l){v^u* ~ u,v*) 

{5^1 + 5^^2){'^u*v* + 2v*v* + u*u* + v*v*) + {5^\ - 5^l){u*u* - «>*) 



B. 



dr ^riTF 



(78) 



17 



Bij = 2 J '^^ \AirF (^^1 + (5<i>2)(2'«iWj + 2viUj + UiUj + ViVj) - ((5$° - (5<i>2)(uiUj - ViVj) 

Starting from Eq. ([77|), one can study both the damping and the frequency shift of the low- lying modes. The 
calculation of the damping rates has been carried out by several authors [pT|-p^. In Refs. |^,|^ the damping of 
the m — and m — 2 mode has been calculated as a function of temperature and found in good agreement with 
experiments QJ^. Concerning the frequency shifts, a calculation based on a method similar to ours has been carried 
out in p7| , but only for quasiclassical modes which satisfy the condition hujfio <^ stf ^ fJ'- In the present work we 
study the frequency shift of the lowest-lying collective modes with exF ~ fi^^ho- These modes have also been studied 
within the dielectric formalism in 24 1. 



B. Frequency shift of the collective modes 

1. Non-resonant contribution 
Similarly to the homogeneous case, the non-resonant contribution to the frequency shift SE is defined as 



<5£;Arfl = 2g2^(l + /f + /j') 



\B,. 



exF 



(79) 



The matrix elements of B and B are given in ([78|). By using the semiclassical approximation ( [70| ) for the quasiparticle 
functions Ui and Vi, one can write SE^ji in the following form 



SE 



NR 



f-TF 



Xo(r) [i^f.f (r, r + s) + k^.^f{v, r + s)] xS(r + s) 



(80) 



where the functions xo(r) are the solutions of (p6[), and in e'**'^['^'(''-'+'^3('')l we have neglected second derivatives of 
the phase (p. The smoothly varying kernels A'f^^^nd ^^j' are defined as 



^2^#(r,r') 



1 + + 2gnTF(r), , ,2gnTF(r') 



6y(r)- 



tTF 



where we have introduced the matrices 



(81) 



fly (r) = 2uj (r)wj (r) -t- 2vi {v)uj (r) 4- {r)uj (r) + Vi {r)vj (r) , 
bij (r) = u, {r)uj (r) - Vi {r)vj (r) . 

In the limit e^F ^ M we can use the following gradient expansion 

Xo(r)i^f,f(r,r + s)xS(r + s)c.i?f^«(r,r) |xo(r)P , 



and 



Xo{r)K^'if{r, r + s)xS(r + s) ^ if-f (r, r) |xo(r)P + ^Xo(r) (s • V)' xS(r) 

Xo(r)(s.V)i^f,f(r,r)(s.V)xS(r) . 



(82) 
(83) 



(84) 



Since the kernel is already zeroth order in exF/fJ'j we can neglect higher order terms in the expansion (^Sf). On 

the contrary, is of order (/x/ctf)^ and we need the expansion (84) to second order in the displacement s. Notice 

also that terms in the expansion (jsj) containing odd powers of s vanish in (|8^) due to geometry. Moreover, we have 
neglected in ( ^4|) second derivatives of the slowly varying functions Ui, u,;. By the replacement and after 

integration by parts, one gets the result 
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5E 



NR 



+ 



pq 



12^2 



yNRi 



(85) 



where k = q + p. In the first term on the r.h.s. of ( pSD the integration over s gives 6qo, while in the second term one 
writes s^g-iqs/n _ _^2y2g-jqs/ri integrates by parts over q. After some algebra one gets the result 



6Enr 



|drgl^(l + 2/«) 



IXo 



2 y^P' |Y7^^ |2"- 6 "tF 



^2„2^2 ,0 



, I |2g'^?^F 1 1 + Vp 



'V ■ 



(86) 



p2 



In the above equation ep(r) and = fp{r), according to (72). We notice that in the homogeneous limit, where 
Yofr) = e'i '"/''/Vy and etf = cbQ, the first two terms on the r.h.s. of (|86| ) coincide with the corresponding terms in 



2. Resonant contribution 



The resonant contribution to Ji? is defined as 

6En = 4g2 - /°) > (87) 

^ ctf + - 

where the matrix elements Aij are given in (^8|) . Following the method used in the analysis of the non- resonant terms 
one has 



SEr _ g 

ETF 2V'^ 



The kernels in the above equation are defined by 
>R 



i?f.(r,r') = -JtJL^ (,^.(r) _ ^I!^.,(r)) (c.,(r') - ^i^d.,(r') 

ETF + - V ~ '^J / V ^3 

K,\M - ^ ^i^rf,(r)^i!^d,(r') , (89) 



— CTF " ■ ■ CTF 

where we have introduced the matrices 



(r) = 2Ui(r)uj(r) + 2w^(r)'yj(r) + w^(r)uj(r) + ?2j(r)'Uj(r) , 

(r) = (r)?2j (r) - Ui {r)vj (r) . (90) 

In the limit ctf IJ- the term in Eq. ( |88| ) containing the kernel K^^j can be treated using the gradient expansion 
(p4). One gets thus 



2^2 



^ 1 drdse-*^-^[^'W-^^«l xo(r)if2''^,(r,r + s)x^(r + s) = 



D J tTF q p ''Iq+Pl 
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In the homogeneous hmit, the first term on the r.h.s. of ( pi| ) coincides with the second term in the square bracket on 
the r.h.s. of (|48|). Moreover, the last term on the r.h.s. of (|8^ ) and ([qi] ) are equal and opposite in sign, and cancel 
out in the sum SEmr + SEn. A similar cancellation is also present in the homogeneous case [see Eqs. (|6|) and 



The contribution to SEf; arising from Kl\j is more delicate. In fact, all terms in the gradient expansion give 
contributions which are of the same order in the limit etf <C /i. Hovever, if we restrict ourselves to modes for which 
V^Xo — const, the expansion (|8^ ) is still appropriate. In fact, higher order terms in ( ^4|) contain derivatives of V^xo 
and vanish for modes with constant Laplacian. To this class of modes belong, for example, all surface modes, for 
which V^xo = 0; E^nd the lowest breathing modes. For the above mentioned term one finds 

I '^^'^^ e— ^[^•W-'^^Wl xo(r)i?fy(r,r + s)xS(r + s) = 

-1 / drdsg^Y.^-^'^'^'^h^K [xo(r)V^xS(r)i?fp,.(r,r) +xo(r)Vxo(r) • Vi?f,,(r,r)] . (92) 
pq 

The Laplacian in momentum space can be easily calculated once the low-q behavior of K^^^ (r, r) has been obtained. 
A straightforward calculation yields: 



4i:'^*.p,-s5;4i:(-|f)';(i 



2 



(93) 



p ^ ^ p 

After some algebra, one obtaines the following result for the contribution to SEb, 

.2 



/ drdse-v[-«-..Ml xo(r)if?.,(r, r + s)xS(r + s) = ^ 



dv\^XQ\''{anTFf'^j dx 



TF 

(w- 1)3/2 ^2u+l 

U + 1 



2 



(94) 



where u — ^1 + t'^{y)x'^ and r(r) = fcsr/gnTF(r)- We stress that the above contribution to the frequency shift is 
peculiar of collective modes with constant Laplacian of harmonically trapped systems in the thermodynamic limit. 
In the homogeneous case, where xo e"i "" and V^xo 7^ const, the contribution of this term is different and is given 
by the first term in the square bracket on the r.h.s. of (Ea). 



3. Results 



We are now in a position to calculate the shift 6E to order g^, by summing the various contributions to the real 
part of Eq. (|77|). One gets the relevant result 



— ^--^^ I dv{anTFf'' (xoV^xS + xSV^Xo) 

+ / dv{a^nTFf'^\xo?G^[T{v)]-^^ f dr {anTF)'^^\^Xo\^GMr)] , 
J muitp J 



(95) 



where ujtf = ^tfI^i and the functions G'i(t), 62 (t) of the local reduced temperature r(r) = / gnxFi^) are 
defined as follows 



Gi(t) 



dx- 



1 



— 1 + u -1 
u u + 1 



(96) 



and 



G2(r) 



dx 



A^/tx- 



V^T^ (u- 1)3/2 /2u+l 



u{u + 1) 



U+1 



^, 1 

dx — 

— 1 u(u + 1) 



(97) 
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where, as usual, u = y/l + t'^x^. The functions G'i(r) and 6*2 (r) are plotted in Fig. 4. Both are positive for any 
value of r. As a consequence, Gi(r) gives an upward shift of the excitation frequency, while G2{t) gives a downward 
shift. The above result holds for collective modes which do not excite the center of mass degrees of freedom. In fact, 
as discussed in Sec. II, the theoretical approach developed in the present work does not describe the center of mass 
motion and, in particular, the dipole mode. 

At r = 0, both Gi and G2 are zero and one is left with the result 



5E 
eTF 



(98) 



which coincides with the findings of Ref. ||2^, obtained from the hydrodynamic theory of superfluids, and of Ref. [ pOf . 
Notice that only non-resonant terms give contribution at T = 0, as a consequence result ( |98|) holds in general for 
collective modes with ctf *C /i, and is not restricted to modes which have constant Laplacian. For the monopole 
(breathing) mode in a spherically symmetric trap (A = 1), characterized by the frequency = V^^ho, one has from 
Eq. (p|) the fractional shift MM 



5ujM 21\/2 



OJM 320C(3) 



(99) 



expressed in terms of the parameter 77. For the m 
frequency given by (|68|), one finds the result [p9 30 



modes in an axially symmetric trap, which have excitation 



320C(3) 



r/V±(A) 



where 



16A2 
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(100) 



(101) 



and the index ± refers to the high (+) and low (— ) m = mode. As discussed in Ref. [g9j, these frequency shifts are 
very small. For 77 = 0.4, which is a typical value for the interaction parameter in experiments, one gets a fractional 
shift of the order of 0.5 %. 

At finite temperature the terms involving the Gi and G2 functions contribute to the frequencyshift, and, differently 
from the T = case, also surface excitations with V^xo = are affected by the correction (|95|). We consider first 
spherically symmetric traps. The monopole oscillation has the form xm oc (r^ — 3i?^/5), where R — yj2^TF{Na) / muj^^ 
is the condensate radius. The temperature dependence of the fractional shift 5lom/^m is given by the equation 



5ujm 21V2 , 
— rj-^ 



UJM 320C(3) ' V N 



1/5 



1 « /"l 

1 + ^ / dxx^/'^VT^(2-5xfGi[T{x)] 
9\/7r Jo 

dxX^^^{l-xf^^G2[Tix)] , 



160 



(102) 



where Nq/N is the equilibrium value of the condensate fraction, which is fixed by the parameter rj and the 
reduced temperature t = T/T^. The argument t{x) of the Gi an d G 2 functions is given by the expression 
t{x) = [{Nq / N)~'^^^t / rj\l / x. In Fig. 6 the monopole frequency shift (102) is shown as a function of the reduced 
temperature t for the value 77 = 0.4 of the interaction parameter. Already at relatively low temperatures, t ~ 0.3, the 
monopole frequency is found to be about 1 % smaller than uom = ^J^^ho- In fact, even for such low temperatures, 
the local reduced temperature is la rge a t the boundary of the condensate, as t{x) ^ 1 if x ^ 0, and t he c ontribution 
of this region dominates the shift (102). If 77 <C t, one can approximate the functions G\ and G2 in (102) with their 
asymptotic behavior for r 3> 1. One has 



Gi(r) 
G2(r) 



Gi(oo)r 
G2(oo)r 



where Gi(oo) — 5^/7r and 
4\/2 



G2(oo) 



3V^ 



dx 



8 + 5a;- 



+ 



V^yr^(l + a;)3 x3/2(l-x2)3/4 



1/4 



(1-^) 

(l + a;)i/4 



22. 



(103) 



(104) 
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In this case the monopole shift is given by 



6uj 



M 



LUM 960C(3) (l-i3)i/5 



[17Gi(oo) - lOGsloo)] 



-1.1 



(l-i3)l/5 ' 



(105) 



where for the condensate fraction we have used the ideal gas law N^/N = 1 — t^. Result ( |105| ) gives a reaso nably 
good approximation to the frequency shift Sujm also when rj ^ t, for example, for r] = 0.4 and t = 0.8 Eq. (|l05| ) gives 
Sujm/^^m — —0.16 and the calculation based on Eq. (102) gives —0.11. 

In a surface mode the oscillation of the condensate has the form xim oc r^Yim{0, 0), and the excitation frequency is 
given by uji ~ y/lmfio. For these modes one finds the following fractional shift 



V2(2? + 3) 
15\AFC(3) ' 



N 



1/5 



-{1 + 1/2) I dxx'^'^{l-xy-^'^G2[T[x)] 



In the limit rj <^t one gets, by using (103), the following result 

V2(2/ + 3) (? + i)r(/ + i) 



30C(3)(; + i) 



IV{1) 



(1-<3)1A 



[2Gi(oo)-G2(cx3)] 



(106) 



(107) 



In the case of axially symmetric traps, the oscillations with m = symmetry have the form Xm=o oc — 2/xtf(s^ — 
2)/[ms^cLi^ + + (s^ — 4)z^], where s — uJm=o/'^± and uJm=o is given by Eq. (pq) for the low and high mode. After 



some algebra one gets the result 



5u!m=0 



21^/2 



Nn 



1/5 



tu^=o 320C(3) ' V N 
16 



v^-^ f±W- 



160 



dxx^/^il-xf/^ G2[Tix)] 



(108) 



9V7r Jo 



4(1 - xf + 3f±{\){7x^ - Ax) 



Gi[t{x) 



where /±(A) is defined in ( |101| ). In the limit i] the above result reduces to 

6.^^o 7^/2^ ,H ^^^^ 3/±(A)]G,(oo) - 10G,(oo)} 



UJm=0 960C(3) (1-^3)1/5 



(109) 



and with excitation energy ujm = y'\m\Lu±, exhibit 



On the contrary, surface excitations of the form Xr> 
the fractional shift (106) with I replaced by \m\. 

In Fig. 5 (Fig. 6) we show the fractional shift of the mode m = low (high) as a function of the reduced temperature 
t = T/T^ and for the value = 0.4 of the interaction parameter. We notice that in the case of the m = high mode 
(Fig. 6), the size of the fractional shift is maximum for spherically symmetric traps (A — 1) and is minimum for 
disk-shaped traps (A ^ 1). On the contrary, for the m = low mode (Fig. 5), jJcj/tjl is maximum for A 1, while it 
is minimum in the A = 1 case. However, by changing the geometry of the trap, the curve of the fractional shift remains 
qualitatively the same, and at intermediate temperatures, T ^ 0.5T^, one finds downward shifts ranging from 1 to 
4% for both the mode to = low and high. In Fig. 7 we show the results for the surface modes to = 2 and to = 4 for 
the same value, rj = 0.4, of the interaction parameter. In the case of surface modes the fractional shift is independent 
of the deformation parameter A and we find that the size of the shift increases by increasing m. An explanation of 
this behavior can be found in the fact that modes with higher to are more localized at the surface of the condensate 
where ksT ^ gno(r), being gno(r) the local chemical potential. Thus, thermal effects are more pronounced for such 
modes. 

Experiments on the temperature dependence of the collective modes have been carried out both at JILA |7j and 
MIT §]. The JILA group has measured, as a function of temperature, the frequency of the to = 2 and to = low 
modes [Q. However, in these experiments, the number of trapped particles is about 10^ and beyond Thomas- Fermi 
effects are expected to play a significant role. Nevertheless, our results for the fractional shift of the m = 2 mode in 
disk-shaped geometries, shown in Fig. 7, both qualitativelly and quantitativelly agree with the observed behavior. 
In the case of the m = low mode other effects, not included in the present analysis, might be responsible for the 
features observed in the experiment. The frequency of collective excitations in the Thomas-Fermi regime has been 
measured by the MIT group for the ?7i = low mode in a cigar-shaped trap M- In Fig- 8 wc show the comparison 
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between the experimental results and our theoretical prediction. The calculation has been carried out with the value 
ry = 0.4 of the interaction parameter, which is close to the experimental conditions of In Fig. 8, the experimental 
data have been plotted as a function of the reduced temperature T/T^ ||6^. This is possible only for temperatures 
above 0.5 /iK, as lower temperatures were not measurable in 



C. Hydrodynamic equations at T = 



At zero temperature, superfluid systems are described by the equations of hydrodynamics (for a general discussion 
see the book ||63| ). These equations involve the total density n of the system and the superfluid velocity Vg, which 
is related to the gradient of the phase of the order parameter. The hydrodynamic picture has been successfully 
employed in |^ to obtain the frequencies of the collective modes in the Thomas-Fermi regime and later in p^] t o 
calculate the corrections to these frequencies due to beyond mean-field effects. We have already verified [see Eq. (|9^)] 
that our perturbation scheme reproduces at T = the results obtained from hydrodynamic theory. However, since 
we start from dynamic equations written in terms of the condensate wavefunction and the noncondensate density, it 
is important to understand whether these equations reduce at zero temperature to the hydrodynamics of superfluids. 

At the level of Gross-Pitaevskii theory the analogy is straightforward . By writing the condensate wavefunction 
in terms of a modulus and a phase <i>(r,t) = Y^no(r, t)e*''''^ one has the following identifications 



5no(r,t) = $o(r)[<5$(r,t)+<5$*(r,t)] , 



(110) 



between the fluctuations of riQ and Lp and the fluctuations of the order parameter. The coupled equations (64), holding 
in the Thomas-Fermi regime, are then equivalent to 



dSriQ 



dt 



TO- 



V • {n.TF^s) = 
f gV(5no = , 



(111) 



where Vj — UVLp/m is the superfluid veloc ity. At T = 0, if one neglects quantum depletion, the condensate density 
coincides with the total density and Eqs. (Ill) coincide with the linearized hydrodynamic equations. The former of 
Eqs. (Ill) corresponds to the equation of continuity and the latter to Euler equation with the pressure P fixed by 
dP/dno = griTF- 

Beyond Gross-Pitaevskii theory one must replace Eqs. by ([76|), which include the corrections to order g^. At 
T = these equations reduce to 



2m 



[5^1 - 5^2) - 2gTO"(5$i - 5^2) 



huj{S^i - 5^2) = "^{griTF + Sn){S^i + 6^2) 

1 -r-^ TO 40 3 



(112) 



2gnTF 



p2 30f 



((5$! -t- 6^2) 



+ gyriTF [(2uiWj + 2viUj + UiUj + ViVj)gij{u)) 

+ (2<w* + 2v*u* + u*u* -f v*v*)g*^{uj) 

where 5^ = ^jl — ^tf{No) is the change in the chemical potential [see Eq. ( |7^)] and the matrices gij{uj) and glj{uj) are 
given in (|7^) with /° = fj = 0. By using the semiclassical approximation (|70|) for the quasiparticle states, the above 
equations are written in terms of the variables ip and hq as 



ihuj Sno{r) = — ?iV • [no(r)vs(r)] — 4gnTF(r)ni-"(i")'5'y9(r) 
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^ X- 



ihiuSno{r + s) 



(113) 



2gnTF{r + s) , 
fly (r + s) H j b^j (r + 



2nTF(r + (r + s)(5(^(r + s) 



and 



ihuj Sip{r) = g ^1 + g-^ XI W 

- 4 y ^ / e--v[..(r)+..(r)l + ,) f (r + s) + 'g"^^(^ + ^) 

^ ei + ej J \ + ej 



(114) 



bij{r + s) 



where the matrices Uij and 6y have been defined in (^2|). By using the gradient expansion employed in the previous 
section for the calculation of the non- resonant contributions to SE, we get the result 



lU! 



(^1 + -^(a^uTFir))'/^^ Snoir) = V • no(r) (^1 + -^{a^TFir))'^'^ v^r) 



and 



ihu/Sip 
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ia'nTFy/^] Snoir) 



(115) 



(116) 



If one takes into account the effect of quantum depletion, the local relation between condensate density and total 
density is given by n = no[l + 8(a^nTF)^''^/(3-\/7r)], and for the fluctuations of the two densities Sn = Sno[l + 
A{a^nTF)^^'^ / Vt^]- Finally, the change in the local chemical potential /i; = gn[l + 32(a^nTF)"^^^/(3\/7r)] induced 
by a density fluctuation is given by the following expression Sndfii/dn = Sngg [1 + 20{a^nTFy^^/yM- It IS now 
straightforward to recognize Eqs. (115), (|116|) as the linearized hydrodynamic equations 



d6n 



V • (nvs) = 



V • 



dn 



Sn] =0 



(117) 



which involve the total density n and the superlfuid velocity Vs 



V. CONCLUDING REMARKS 



In this paper we have studied the coUisionless collective modes of a dilute Bose gas beyond the Gross-Pitaevskii 
theory. In particular, for harmonically trapped systems in the thermodynamic limit, we have calculated the corrections 
to the excitation frequencies of the low-lying collective modes. We find that, not far below the Bose-Einstein transition 
temperature, the fractional frequency shift is of the order of few percent for typical experimental conditions and can 
be measured. A direct comparison with experimental data obtained by the MIT group with large condensates looks 
very good. Similarly to what happened to Gross-Pitaevskii theory, the study of collective excitations can become a 
useful bench-mark also for theories beyond mean-field approximation. 
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FIG. 2. Dimensionless function _F as a function of tire reduced temperature r — ksT/gno (solid line). The asymptotic 
behaviors for r ^ 1 (dashed line) and for t ^ 1 (long-dashed line) are also reported. 



FIG. 5. Fractional shift of the m — Q low mode as a function of the reduced temperature T/T^. The value of the interaction 
parameter is ?7 = 0.4. For a spherical trap (A = 1) this mode corresponds to the I = 2 quadrupole mode. 

FIG. 6. Fractional shift of the m = high mode as a function of the reduced temperature T/T^. The value of the interaction 
parameter is r; = 0.4. For a spherical trap (A = 1) this mode corresponds to the monopole (breathing) mode. 

FIG. 7. Fractional shift of the m — 2 and m — A surface modes as a function of the reduced temperature T/T°. The value 
of the interaction parameter is rj = 0.4. 

FIG. 8. Temperature dependence of the frequency of the m = low mode. The experimental points are taken from The 
theoretical calculation (solid line) corresponds to cigar-shaped traps and to the value ri = 0.4 of the interaction parameter. The 
dashed line corresponds to the hydrodynamic prediction u = ^hjlvz from Eq. (S). 




FIG. 1. Dimensionless function _ff as a function of the reduced temperature t = ksT /gno. 



FIG. 3. Dimensionless function G as a function of the reduced temperature r = fesT/gno- 



FIG. 4. Dimensionless functions Gl and G2 as a function of the reduced temperature r = kBT/gno. 




26 



Fig. 1 




X 



-2 



0.0 



0.5 



1.0 



1.5 



2.0 



2.5 



3.0 



Fig. 4 




X 



Fig. 5 




Fig. 6 




Fig. 7 




Fig. 8 



N 



28 



27 - 



26 - 



25 

(D 

1 24 



23 



22 




T/T 



